#include "easyrun.hpp"
#include "option_parse.hpp"
#include "call_rates.hpp"
#include "snp_call.hpp"
#include "consensus.hpp"

void do_easyrun(int argc, char **argv)
{
    prog_options options = parse_easyrun_options(argc, argv);

    /* load the lengths of all the genomes. */
    map_str_ui_t ref_seq_len = build_fasta_index(options.fasta_file);

    unsigned int read_length = get_read_length(options.map_file.c_str());
    mapping_type t = get_mapping_type(options.map_file.c_str());

    if(t & COLOR )
    {
        // TODO: UGLY. Find a better way to do this.
        options.old_bin_cnt = 36;

        map_str_vafl_t subs, totl;
        cand_blocks_t cand_blocks;
        cand_vector_t snp_cands;

        find_subs_totl(options, subs, totl, ref_seq_len);
        find_snp_cands(options, cand_blocks, snp_cands, subs, totl);
        
        std::map<std::string, std::string> seqs;
        parse_fasta(options.fasta_file.c_str(), seqs);
        attach_seqs_to_snp_cands(seqs, snp_cands);

        array_type5 call_counts, call_rates;
        
        find_call_rates(options, call_counts, cand_blocks, snp_cands);
        boost::multi_array<int, 2> qual_bins;

        // rebin quals.
        rebin_quals(call_counts, read_length, options.old_bin_cnt, options.new_bin_cnt, qual_bins);
        
        // use the new bin boundaries to collapse the call_rates matrix.
        rebin_calls(call_counts, qual_bins, read_length, options.old_bin_cnt, options.new_bin_cnt, call_rates);

        // convert the counts to probabilities.
        convert_call_rates_to_probs(call_rates, read_length, options.new_bin_cnt);

        color_call_vector_t snp_calls;
        for(unsigned int i = 0; i < snp_cands.size(); i++)
        {
            SnpRecord &snp = snp_cands[i];
            snp_calls.push_back(ColorSnpCall(snp.id, snp.chr, snp.start, snp.end, snp.seq, snp.subs, snp.totl));
        }
        
        call_snps(options, cand_blocks, snp_calls, call_rates, read_length, qual_bins);

        multiply_col_snp_calls_by_priors(options, snp_calls);

        if(options.cons_file != "")
        {
            std::map<std::string, std::string> genome_seqs;
            parse_fasta(options.fasta_file.c_str(), genome_seqs);

            change_col_genome(options, snp_calls, genome_seqs);

            std::ofstream ofs(options.cons_file.c_str());

            write_fasta(ofs, genome_seqs);

            ofs.close();
        }

        print_col_snps(options, snp_calls);
    } 
    
    else // SEQUENCE SPACE
    {
        // TODO: UGLY. Find a better way to do this.
        options.old_bin_cnt = 96;

        map_str_vafl_t subs, totl;
        cand_blocks_t cand_blocks;
        cand_vector_t snp_cands;

        find_seq_subs_totl(options, subs, totl, ref_seq_len);
        find_seq_snp_cands(options, cand_blocks, snp_cands, subs, totl);
        
        std::map<std::string, std::string> seqs;
        parse_fasta(options.fasta_file.c_str(), seqs);
        attach_seqs_to_snp_cands(seqs, snp_cands);

        array_type4 call_counts, call_rates;
        
        find_seq_call_rates(options, call_counts, cand_blocks, snp_cands);
        boost::multi_array<int, 2> qual_bins;
        
        // rebin quals.
        rebin_seq_quals(call_counts, read_length, options.old_bin_cnt, options.new_bin_cnt, qual_bins);
        
        // use the new bin boundaries to collapse the call_rates matrix.
        rebin_seq_calls(call_counts, qual_bins, read_length, options.old_bin_cnt, options.new_bin_cnt, call_rates);

        // convert the counts to probabilities.
        convert_seq_call_rates_to_probs(call_rates, read_length, options.new_bin_cnt);

        seq_call_vector_t snp_calls;
        for(unsigned int i = 0; i < snp_cands.size(); i++)
        {
            SnpRecord &snp = snp_cands[i];
            snp_calls.push_back(SeqSnpCall(snp.id, snp.chr, snp.start, snp.end, snp.seq, snp.subs, snp.totl));
        }
        call_seq_snps(options, cand_blocks, snp_calls, call_rates, read_length, qual_bins);

        multiply_seq_snp_calls_by_priors(options, snp_calls);
        
        if(options.cons_file != "")
        {
            std::map<std::string, std::string> genome_seqs;
            parse_fasta(options.fasta_file.c_str(), genome_seqs);
            
            for(unsigned int i = 1; i < snp_calls.size(); i++)
            {
                SeqSnpCall &snp = snp_calls[i];

                if( snp.seq.empty() )
                    continue;

                unsigned int max_idx = std::max_element(snp.probs.begin(), snp.probs.end()) - snp.probs.begin();
                genome_seqs[snp.chr][snp.start] = snp.blocks[max_idx].base;
            }

            std::ofstream ofs(options.cons_file.c_str());
            write_fasta(ofs, genome_seqs);
            ofs.close();
        }

        print_seq_snps(options, snp_calls);
    }

}
